Mean Sea Surface Temperature#

../../../_images/d5e2a63926bbfe0c9d8d451cea54d8ea32aec49f38d4a415f62a1995c8580e83.png
../../../_images/eb16e9afefb2e55321589a5dd93c51ac39081c768d8c6a75cf0860ab3d67c030.png

Figure. Sea Surface Temperature (SST) Change from satellite. The map (top) shows the change in mean SST (°C) in the vicinity of Palau over the period 1982–2024 from the NOAA OISSTv2 satellite. The grey line is the Palau EEZ. The bar plot (bottom) shows the mean SST averaged over the area within the top plot. The trend is statistically significant (p < 0.05). The colored dots represent the 10 warmest years on record.

Setup#

First, we need to import all the necessary libraries. Some of them are specifically developed to handle the download and plotting of the data and are hosted at the indicators set-up repository in GitHub

https://www.ncei.noaa.gov/products/optimum-interpolation-sst

Hide code cell source
import warnings
warnings.filterwarnings("ignore")

import sys
import os
import os.path as op
import xarray as xr
import geopandas as gpd
import pandas as pd
import cartopy.crs as ccrs
import numpy as np

import matplotlib.pyplot as plt
from myst_nb import glue 

sys.path.append("../../../../indicators_setup")
from ind_setup.plotting import plot_base_map, plot_bar_probs, plot_map_subplots, fontsize
from ind_setup.plotting_int import plot_timeseries_interactive, plot_oni_index_th

sys.path.append("../../../functions")
from data_downloaders import download_oni_index
lon_site, lat_site = 134.368203,7.322074

#Area of interest
lon_range  = [129.4088, 137.0541]
lat_range = [1.5214, 11.6587]
shp_f = op.join(os.getcwd(), '..', '..','..', 'data/Palau_EEZ/pw_eez_pol_april2022.shp')
shp_eez = gpd.read_file(shp_f)
path_data = "../../../data"
path_figs = "../../../matrix_cc/figures"
data = xr.load_dataset(op.join(path_data, 'sst_daily_1981_2024_palau.nc'))
dataset_id = 'sst'

Analysis#

fig, ax = plot_base_map(shp_eez = shp_eez, figsize = [10, 6])
im = ax.pcolor(data.lon, data.lat, data.mean(dim='time')[dataset_id], transform=ccrs.PlateCarree(), 
                cmap = 'hot_r', vmin = np.percentile(data.mean(dim = 'time')[dataset_id], 1), 
                vmax = np.percentile(data.mean(dim = 'time')[dataset_id], 99))
ax.set_extent([lon_range[0], lon_range[1], lat_range[0], lat_range[1]], crs=ccrs.PlateCarree())
plt.colorbar(im, ax=ax, label='SST (ºC)')

glue("average_map", fig, display=False)
plt.savefig(op.join(path_figs, 'F12_SST_map.png'), dpi=300, bbox_inches='tight')
../../../_images/d5e2a63926bbfe0c9d8d451cea54d8ea32aec49f38d4a415f62a1995c8580e83.png

Annual average#

data_y = data.resample(time='1YE').mean()
im = plot_map_subplots(data_y, dataset_id, shp_eez = shp_eez, cmap = 'hot_r', 
                  vmin = np.nanpercentile(data_y.min(dim = 'time')[dataset_id], 1), 
                  vmax = np.nanpercentile(data_y.max(dim = 'time')[dataset_id], 99),
                  cbar = 1, cbar_pad = 0.05)
../../../_images/8ecca7714035b3c56d8a3efa9234cad327a296a0ca03e84db8c44732652f0011.png

Annual anomaly#

data_an = data_y - data.mean(dim='time')
fig = plot_map_subplots(data_an, dataset_id, shp_eez = shp_eez, cmap='RdBu_r',
                        vmin=-2, vmax=2, cbar = 1, cbar_pad = 0.05)
../../../_images/865ee4492f852129e25cbf141ca2dcd03cb799e2d92494d9390535563bc0a6c4.png

Average Area#

data_mean = data.mean(dim = ['lon', 'lat']).to_dataframe()
datag = data_mean.groupby(data_mean.index.year).max()
datag.index = pd.to_datetime(datag.index, format = '%Y')
datag['sst_ref'] = datag['sst'] - datag.loc['1961':'1990'].sst.mean()
nevents = 10
data_mean = data.mean(dim = ['lon', 'lat']).to_dataframe()
data_mean = data_mean.resample('YE').mean()
top_10 = data_mean.sort_values(by='sst', ascending=False).head(nevents)


fig, ax, trend = plot_bar_probs(x=data_mean.index.year, y=data_mean.sst, trendline=True,
                                y_label='SST [°C]', figsize=[15, 4], return_trend=True)
ax.set_ylim(data_mean.sst.min()-.2, data_mean.sst.max()+.2)
im = ax.scatter(top_10.index.year, top_10.sst, 
                c=top_10.sst.values, s=100, ec = 'pink', cmap='rainbow', label='Top 10 warmest years')
plt.colorbar(im, pad = 0).set_label('Absolute SST [°C]', fontsize=fontsize)

ax.set_title('Annual Mean SST', fontsize=15);
glue("average_timeseries", fig, display=False)
plt.savefig(op.join(path_figs, 'F12_SST_trends_Annomalies.png'), dpi=300, bbox_inches='tight')
../../../_images/eb16e9afefb2e55321589a5dd93c51ac39081c768d8c6a75cf0860ab3d67c030.png

Given point#

loc = [7.37, 134.7]
dict_plot = [{'data' : data.sel(lon=loc[1], lat=loc[0], method='nearest').to_dataframe(), 
              'var' : dataset_id, 'ax' : 1, 'label' : 'SST [ºC]'},]
fig, ax = plot_base_map(shp_eez = shp_eez, figsize = [10, 6])
ax.set_extent([lon_range[0], lon_range[1], lat_range[0], lat_range[1]], crs=ccrs.PlateCarree())
ax.plot(loc[1], loc[0], '*', markersize = 12, color = 'royalblue', transform=ccrs.PlateCarree(), label = 'Location Analysis')
ax.legend()
<matplotlib.legend.Legend at 0x181b4e510>
../../../_images/038a2757803511e8b284fa11a57d0ba084269dbabe56e09215299e69572e231d.png
fig = plot_timeseries_interactive(dict_plot, trendline=True, scatter_dict = None, figsize = (25, 12), label_yaxes = 'SST [ºC]');

ONI index analysis#

p_data = 'https://psl.noaa.gov/data/correlation/oni.data'
df1 = download_oni_index(p_data)
lims = [-.5, .5]
plot_oni_index_th(df1, lims = lims)

Group by ONI category

data_xr = xr.load_dataset(op.join(path_data, 'sst_daily_1981_2024_palau.nc')).rename({'lat':'latitude', 'lon':'longitude'}).resample(time='1M').mean()
data_xr['time'] = data_xr.time.values.astype('datetime64[M]')

data_xr['ONI'] = (('time'), df1.iloc[np.intersect1d(data_xr.time, df1.index, return_indices=True)[2]].ONI.values)
data_xr['ONI_cat'] = (('time'), np.where(data_xr.ONI < lims[0], -1, np.where(data_xr.ONI > lims[1], 1, 0)))
data_oni = data_xr.groupby('ONI_cat').mean()
data_oni['labels'] = (('ONI_cat'),['La Niña', 'Neutral', 'El Niño'])
dataset_id = 'sst'
fig = plot_map_subplots(data_oni, dataset_id, shp_eez = shp_eez, cmap = 'hot_r', 
                  vmin = np.percentile(data_xr.mean(dim = 'time')[dataset_id], 0) - .25, 
                  vmax = np.percentile(data_xr.mean(dim = 'time')[dataset_id], 100) + .25,
                  sub_plot= [1, 3], figsize = (15, 8), cbar = True, cbar_pad = .1)

plt.savefig(op.join(path_figs, 'F14_SST_ENSO.png'), dpi=300, bbox_inches='tight')
../../../_images/e3f56ecbd4b2cada7c4a664bd4c1b6e03ee614733842c3ba96ad845abf3c67e2.png
data_an = data_oni - data_xr.mean(dim='time')
data_an['labels'] = (('ONI_cat'),['La Niña', 'Neutral', 'El Niño'])
fig = plot_map_subplots(data_an, dataset_id, shp_eez = shp_eez, cmap='RdBu_r', vmin=-.4, vmax=.4,
                  sub_plot= [1, 3], figsize = (15, 8), cbar = True, cbar_pad = 0.1)
../../../_images/6ae2c6bcb79a79ce48d7b6e195bef095ae22793262242eabed64c85077b97424.png